module radar_bin_mod

  use container
  use flogger
  use data_translators_utils_mod
  use radar_mod
  use string

  implicit none

  private

  public radar_bin_read

  integer, parameter :: max_radial_size = 370

contains

  subroutine radar_bin_read(file_path, stations, records)

    character(*), intent(in) :: file_path
    type(hash_table_type), intent(inout) :: stations
    type(linked_list_type), intent(inout) :: records

    ! Common header fields
    integer(4) magic_number                       ! 01
    integer(2) major_version                      ! 02
    integer(2) minor_version                      ! 03
    integer(4) generic_type                       ! 04
    integer(4) product_type                       ! 05
    character(16) reserved1                       ! 06
    ! Station fields
    character(8) site_id                          ! 01
    character(32) site_name                       ! 02
    real(4) site_lat                              ! 03
    real(4) site_lon                              ! 04
    integer(4) antenna_hgt                        ! 05
    integer(4) site_alt                           ! 06
    real(4) freq                                  ! 07
    real(4) beam_width_h                          ! 08
    real(4) beam_width_v                          ! 09
    integer(4) rda_version                        ! 10
    integer(2) radar_type                         ! 11
    character(54) reserved2                       ! 12
    ! Task fields
    character(32) task_name                       ! 01
    character(128) task_desc                      ! 02
    integer(4) polar_type                         ! 03
    integer(4) scan_type                          ! 04
    integer(4) pulse_width                        ! 05
    integer(4) scan_start_time                    ! 06
    integer(4) cut_number                         ! 07
    real(4) noise_h                               ! 08
    real(4) noise_v                               ! 09
    real(4) calib_h                               ! 10
    real(4) calib_v                               ! 11
    real(4) noise_temp_h                          ! 12
    real(4) noise_temp_v                          ! 13
    real(4) zdr_calib                             ! 14
    real(4) phidp_calib                           ! 15
    real(4) ldr_calib                             ! 16
    character(40) reserved3                       ! 17
    ! Scan fields
    integer(4) process_mode                       ! 01
    integer(4) wave_form                          ! 02
    real(4) prf1                                  ! 03
    real(4) prf2                                  ! 04
    integer(4) dealias_mode                       ! 05
    real(4) azimuth                               ! 06
    real(4) elevation                             ! 07
    real(4) start_angle                           ! 08
    real(4) end_angle                             ! 09
    real(4) angular_res                           ! 10
    real(4) scan_speed                            ! 11
    integer(4) log_res                            ! 12
    integer(4) doppler_res                        ! 13
    integer(4) max_range_1                        ! 14
    integer(4) max_range_2                        ! 15
    integer(4) start_range                        ! 16
    integer(4) sample_1                           ! 17
    integer(4) sample_2                           ! 18
    integer(4) phase_mode                         ! 19
    real(4) atm_loss                              ! 20
    real(4) nyq_speed                             ! 21
    integer(8) moments_mask                       ! 22
    integer(8) moments_size_mask                  ! 23
    integer(4) misc_filter_mask                   ! 24
    real(4) sqi_thresh                            ! 25
    real(4) sig_thresh                            ! 26
    real(4) csr_thresh                            ! 27
    real(4) log_thresh                            ! 28
    real(4) cpa_thresh                            ! 29
    real(4) pmi_thresh                            ! 30
    real(4) dplog_thresh                          ! 31
    character(4) reserved4                        ! 32
    integer(4) dbt_mask                           ! 33
    integer(4) dbz_mask                           ! 34
    integer(4) vel_mask                           ! 35
    integer(4) spec_width_mask                    ! 36
    integer(4) dp_mask                            ! 37
    character(12) reserved5                       ! 38
    character(4) reserved6                        ! 39
    integer(4) scan_dir                           ! 40
    integer(2) gnd_clut_classi_type               ! 41
    integer(2) gnd_clut_filter_type               ! 42
    integer(2) gnd_clut_filter_notch_width        ! 43
    integer(2) gnd_clut_filter_window             ! 44
    character(72) reserved7                       ! 45
    ! Radial header fields
    integer(4) radial_state                       ! 01
    integer(4) spot_blank                         ! 02
    integer(4) sequence_number                    ! 03
    integer(4) radial_number                      ! 04
    integer(4) elevation_number                   ! 05
    ! real(4) azimuth                             ! 06
    ! real(4) elevation                           ! 07
    integer(4) seconds                            ! 08
    integer(4) microseconds                       ! 09
    integer(4) data_length                        ! 10
    integer(4) moment_number                      ! 11
    character(2) reserved8                        ! 12
    integer(2) est_noise_h                        ! 13
    integer(2) est_noise_v                        ! 14
    character(14) reserved9                       ! 15
    ! Radial data header fields
    integer(4) data_type                          ! 01
    integer(4) scale                              ! 02
    integer(4) offset                             ! 03
    integer(2) bin_length                         ! 04
    integer(2) flags                              ! 05
    integer(4) length                             ! 06
    character(12) reserved10                      ! 07
    ! Radial data fields
    integer(1), allocatable :: bins_1(:)
    integer(2), allocatable :: bins_2(:)

    type(datetime_type) :: time
    type(radar_station_type), pointer :: station
    type(radar_record_type) record
    integer iscan, imoment, iradial, ibin
    character(8) name

    call log_notice('Reading ' // trim(file_path) // '...')

    open(10, file=file_path, status='old', form='unformatted', access='stream')

    ! Common fields
    read(10) magic_number                         ! 01
    read(10) major_version                        ! 02
    read(10) minor_version                        ! 03
    read(10) generic_type                         ! 04
    read(10) product_type                         ! 05
    read(10) reserved1                            ! 06
    if (magic_number /= 1297371986) then
      write(*, *) '[Error]: Invalid data file!'
      stop 1
    end if
    if (generic_type /= 1) then
      write(*, *) '[Error]: Invalid file type!'
      write(*, *) 'generic_type = ', generic_type
      stop 1
    end if
    ! Station fields
    read(10) site_id                              ! 01
    read(10) site_name                            ! 02
    read(10) site_lat                             ! 03
    read(10) site_lon                             ! 04
    read(10) antenna_hgt                          ! 05
    read(10) site_alt                             ! 06
    read(10) freq                                 ! 07
    read(10) beam_width_h                         ! 08
    read(10) beam_width_v                         ! 09
    read(10) rda_version                          ! 10
    read(10) radar_type                           ! 11
    read(10) reserved2                            ! 12

    if (.not. stations%hashed(site_id)) then
      allocate(station)
      call station%init(site_id, site_name, site_lon, site_lat, real(site_alt), real(antenna_hgt), freq)
      call stations%insert(site_id, station)
      deallocate(station)
    end if
    select type (value => stations%value(site_id))
    type is (radar_station_type)
      station => value
    end select

    record%seq_id = records%size
    record%station => station

    ! Task fields
    read(10) task_name                            ! 01
    read(10) task_desc                            ! 02
    read(10) polar_type                           ! 03
    read(10) scan_type                            ! 04
    read(10) pulse_width                          ! 05
    read(10) scan_start_time                      ! 06
    read(10) cut_number                           ! 07
    read(10) noise_h                              ! 08
    read(10) noise_v                              ! 09
    read(10) calib_h                              ! 10
    read(10) calib_v                              ! 11
    read(10) noise_temp_h                         ! 12
    read(10) noise_temp_v                         ! 13
    read(10) zdr_calib                            ! 14
    read(10) phidp_calib                          ! 15
    read(10) ldr_calib                            ! 16
    read(10) reserved3                            ! 17
    ! Check inputs.
    call check_range_with_error(polar_type, 1, 4, 'polar_type', __FILE__, __LINE__)
    call check_range_with_error(scan_type, 0, 6, 'scan_type', __FILE__, __LINE__)
    call check_range_with_error(pulse_width, 1, 1000000, 'pulse_width', __FILE__, __LINE__)
    call check_range_with_error(cut_number, 1, 256, 'cut_number', __FILE__, __LINE__)
    call check_range_with_error(noise_h, -100.0, 0.0, 'noise_h', __FILE__, __LINE__)
    call check_range_with_error(noise_v, -100.0, 0.0, 'noise_v', __FILE__, __LINE__)
    call check_range_with_error(calib_h, 0.0, 200.0, 'calib_h', __FILE__, __LINE__)
    call check_range_with_warning(calib_v, 0.0, 200.0, 'calib_v', __FILE__, __LINE__)
    call check_range_with_error(noise_temp_h, 0.0, 800.0, 'noise_temp_h', __FILE__, __LINE__)
    call check_range_with_error(noise_temp_v, 0.0, 800.0, 'noise_temp_v', __FILE__, __LINE__)
    call check_range_with_error(zdr_calib, -10.0, 10.0, 'zdr_calib', __FILE__, __LINE__)
    call check_range_with_warning(phidp_calib, -180.0, 180.0, 'phidp_calib', __FILE__, __LINE__)
    call check_range_with_error(ldr_calib, -60.0, 0.0, 'ldr_calib', __FILE__, __LINE__)
    ! Initialize record object.
    call time%init(timestamp=scan_start_time)
    call record%init(task_name, time, cut_number)
    select case (scan_type)
    case (0)
      record%scan_type = 'PPI'
    case default
      call log_error('Radar scan_type is ' // to_str(scan_type) // '!')
    end select
    ! Scan fields
    do iscan = 1, cut_number
      read(10) process_mode                       ! 01
      read(10) wave_form                          ! 02
      read(10) prf1                               ! 03
      read(10) prf2                               ! 04
      read(10) dealias_mode                       ! 05
      read(10) azimuth                            ! 06
      read(10) elevation                          ! 07
      read(10) start_angle                        ! 08
      read(10) end_angle                          ! 09
      read(10) angular_res                        ! 10
      read(10) scan_speed                         ! 11
      read(10) log_res                            ! 12
      read(10) doppler_res                        ! 13
      read(10) max_range_1                        ! 14
      read(10) max_range_2                        ! 15
      read(10) start_range                        ! 16
      read(10) sample_1                           ! 17
      read(10) sample_2                           ! 18
      read(10) phase_mode                         ! 19
      read(10) atm_loss                           ! 20
      read(10) nyq_speed                          ! 21
      read(10) moments_mask                       ! 22
      read(10) moments_size_mask                  ! 23
      read(10) misc_filter_mask                   ! 24
      read(10) sqi_thresh                         ! 25
      read(10) sig_thresh                         ! 26
      read(10) csr_thresh                         ! 27
      read(10) log_thresh                         ! 28
      read(10) cpa_thresh                         ! 29
      read(10) pmi_thresh                         ! 30
      read(10) dplog_thresh                       ! 31
      read(10) reserved4                          ! 32
      read(10) dbt_mask                           ! 33
      read(10) dbz_mask                           ! 34
      read(10) vel_mask                           ! 35
      read(10) spec_width_mask                    ! 36
      read(10) dp_mask                            ! 37
      read(10) reserved5                          ! 38
      read(10) reserved6                          ! 39
      read(10) scan_dir                           ! 40
      read(10) gnd_clut_classi_type               ! 41
      read(10) gnd_clut_filter_type               ! 42
      read(10) gnd_clut_filter_notch_width        ! 43
      read(10) gnd_clut_filter_window             ! 44
      read(10) reserved7                          ! 45

      call check_range_with_error(process_mode, 1, 2, 'process_mode', __FILE__, __LINE__)
      call check_range_with_warning(wave_form, 0, 6, 'wave_form', __FILE__, __LINE__)
      call check_range_with_error(prf1, 1.0, 10000.0, 'prf1', __FILE__, __LINE__)
      if (dealias_mode > 1) call check_range_with_error(prf2, 1.0, 10000.0, 'prf2', __FILE__, __LINE__)
      call check_range_with_error(dealias_mode, 1, 4, 'dealias_mode', __FILE__, __LINE__)
      if (process_mode == 1) then
        call check_range_with_error(angular_res, 0.0, 2.0, 'angular_res', __FILE__, __LINE__)
      end if
      call check_range_with_error(scan_speed, 0.0, 100.0, 'scan_speed', __FILE__, __LINE__)
      call check_range_with_error(log_res, 0, 5000, 'log_res', __FILE__, __LINE__)
      call check_range_with_error(doppler_res, 0, 5000, 'doppler_res', __FILE__, __LINE__)
      call check_range_with_error(max_range_1, 0, 500000, 'max_range_1', __FILE__, __LINE__)
      if (dealias_mode > 1) call check_range_with_error(max_range_2, 0, 500000, 'max_range_2', __FILE__, __LINE__)
      call check_range_with_error(start_range, 0, 500000, 'start_range', __FILE__, __LINE__)
      call check_range_with_error(sample_1, 2, 512, 'sample_1', __FILE__, __LINE__)
      if (dealias_mode > 1) call check_range_with_error(sample_2, 2, 512, 'sample_2', __FILE__, __LINE__)
      call check_range_with_warning(phase_mode, 1, 3, 'phase_mode', __FILE__, __LINE__)
      call check_range_with_error(atm_loss, 0.0, 10.0, 'atm_loss', __FILE__, __LINE__)
      call check_range_with_error(nyq_speed, 0.0, 100.0, 'nyq_speed', __FILE__, __LINE__)
      call check_range_with_error(sqi_thresh, 0.0, 1.0, 'sqi_thresh', __FILE__, __LINE__)
      call check_range_with_error(sig_thresh, 0.0, 20.0, 'sig_thresh', __FILE__, __LINE__)
      call check_range_with_error(csr_thresh, 0.0, 100.0, 'csr_thresh', __FILE__, __LINE__)
      call check_range_with_error(log_thresh, 0.0, 20.0, 'log_thresh', __FILE__, __LINE__)
      call check_range_with_error(cpa_thresh, 0.0, 100.0, 'cpa_thresh', __FILE__, __LINE__)
      call check_range_with_error(pmi_thresh, 0.0, 1.0, 'pmi_thresh', __FILE__, __LINE__)
      call check_range_with_error(dplog_thresh, 0.0, 100.0, 'dplog_thresh', __FILE__, __LINE__)
      call check_range_with_error(dbt_mask, 0, 1, 'dbt_mask', __FILE__, __LINE__)
      call check_range_with_error(dbz_mask, 0, 1, 'dbz_mask', __FILE__, __LINE__)
      call check_range_with_error(vel_mask, 0, 1, 'vel_mask', __FILE__, __LINE__)
      call check_range_with_error(spec_width_mask, 0, 1, 'spec_width_mask', __FILE__, __LINE__)
      call check_range_with_warning(dp_mask, 0, 1, 'dp_mask', __FILE__, __LINE__)
      call check_range_with_error(scan_dir, 0, 1, 'scan_dir', __FILE__, __LINE__)
      call check_range_with_error(gnd_clut_classi_type, 1_2, 4_2, 'gnd_clut_classi_type', __FILE__, __LINE__)
      call check_range_with_error(gnd_clut_filter_type, 0_2, 5_2, 'gnd_clut_filter_type', __FILE__, __LINE__)
      call check_range_with_error(gnd_clut_filter_notch_width, 0_2, 10_2, 'gnd_clut_filter_notch_width', __FILE__, __LINE__)
      call check_range_with_error(gnd_clut_filter_window, 0_2, 4_2, 'gnd_clut_filter_window', __FILE__, __LINE__)
      call record%scans(iscan)%init(elevation, real(log_res), real(doppler_res), nyq_speed, max_radial_size)
    end do

    do iscan = 1, size(record%scans)
      do while (.true.)
        ! Radial header fields
        read(10) radial_state                     ! 01
        read(10) spot_blank                       ! 02
        read(10) sequence_number                  ! 03
        read(10) radial_number                    ! 04
        read(10) elevation_number                 ! 05
        read(10) azimuth                          ! 06
        read(10) elevation                        ! 07
        read(10) seconds                          ! 08
        read(10) microseconds                     ! 09
        read(10) data_length                      ! 10
        read(10) moment_number                    ! 11
        read(10) reserved8                        ! 12
        read(10) est_noise_h                      ! 13
        read(10) est_noise_v                      ! 14
        read(10) reserved9                        ! 15
        call check_range_with_error(radial_state, 0, 6, 'radial_state', __FILE__, __LINE__)
        call check_range_with_error(spot_blank, 0, 1, 'spot_blank', __FILE__, __LINE__)
        call check_range_with_error(azimuth, 0.0, 360.0, 'azimuth', __FILE__, __LINE__)
        call check_range_with_error(elevation, -2.0, 90.0, 'elevation', __FILE__, __LINE__)
        call check_range_with_error(moment_number, 1, 64, 'moment_number', __FILE__, __LINE__)

        if (radial_state == 0 .or. radial_state == 3) then
          iradial = 1
        end if

        if (iradial > max_radial_size) then
          call log_error('Insufficient radials array size, enlarge max_radial_size!', __FILE__, __LINE__)
        end if

        ! print *, '==> azimuth ', azimuth, ' elevation ', elevation, seconds+microseconds*1.0d-6
        call time%init(timestamp=seconds+microseconds*1.0d-6)
        call record%scans(iscan)%radials(iradial)%init(time, azimuth, elevation, moment_number)
        do imoment = 1, size(record%scans(iscan)%radials(iradial)%moments)
          read(10) data_type                      ! 01
          read(10) scale                          ! 02
          read(10) offset                         ! 03
          read(10) bin_length                     ! 04
          read(10) flags                          ! 05
          read(10) length                         ! 06
          read(10) reserved10                     ! 07
          call check_range_with_error(data_type, 1, 64, 'data_type', __FILE__, __LINE__)
          select case (data_type)
          case (1)
            name = 'dBT'
          case (2)
            name = 'dBZ'
          case (3)
            name = 'V'
          case (4)
            name = 'W'
          case (5)
            name = 'SQI'
          case (6)
            name = 'CPA'
          case (7)
            name = 'ZDR'
          case (8)
            name = 'LDR'
          case (9)
            name = 'CC'
          case (10)
            name = 'DP'
          case (11)
            name = 'KDP'
          case (12)
            name = 'CP'
          case (14)
            name = 'HCL'
          case (15)
            name = 'CF'
          case (16)
            name = 'SNRH'
          case (17)
            name = 'SNRV'
          case (19)
            name = 'POTS'
          case (21)
            name = 'COP'
          case (26)
            name = 'VELSZ'
          case (27)
            name = 'DR'
          case (32)
            name = 'Zc'
          case (33)
            name = 'Vc'
          case (34)
            name = 'Wc'
          case (35)
            name = 'ZDRc'
          case default
            name = 'reserved'
          end select
          length = length / bin_length
          call check_range_with_error(bin_length, 1_2, 2_2, 'bin_length', __FILE__, __LINE__)

          call record%scans(iscan)%radials(iradial)%moments(imoment)%init(name, real(scale), real(offset), length)

          select case (bin_length)
          case (1)
            allocate(bins_1(length))
            read(10) bins_1
            do ibin = 1, length
              record%scans(iscan)%radials(iradial)%moments(imoment)%bins(ibin) = bins_1(ibin) * scale + offset
            end do
            deallocate(bins_1)
          case (2)
            allocate(bins_2(length))
            read(10) bins_2
            do ibin = 1, length
              record%scans(iscan)%radials(iradial)%moments(imoment)%bins(ibin) = bins_2(ibin) * scale + offset
            end do
            deallocate(bins_2)
          end select
        end do
        if (radial_state == 2 .or. radial_state == 4) then
          record%scans(iscan)%nradial = iradial
          exit
        end if
        iradial = iradial + 1
      end do
      call log_print('Elevationn ' // to_str(record%scans(iscan)%elevation, 3, 7) // ' with ' // to_str(moment_number) // ' data types.')
    end do

    call records%insert(record)

    close(10)

  end subroutine radar_bin_read

end module radar_bin_mod
